% 12/19/2023:
% this version is written for the single channel (IOS) opto stim experiment

%% load in previously selected referencing information:
[fn,fd] = uigetfile('*.mat','select referencing information .mat file');
cd(fd);
load(fn);

%% load in vessel mask:
load(vesselmaskname,'vesselMasks','avios','refim_windowmask','micronsPerPixel')     %avios is what we're going to be registering everything to; there is ONE of these per mouse 
refim = avios;              %need to make a reference ios image for each session

%% load in background:
cd(bgdirectory);
bg_ios = loadRawStack_v1('ackground',localimagesize);

load(bgmetaname,'numrepeats','numframes')

bg_ios = reshape(bg_ios,size(bg_ios,1),size(bg_ios,2),numframes,numrepeats);

% bguid = [reshape(rectPositions,2,size(rectPositions,3));movindicator];
% [uid_bg,~,bgID_idx] = unique(bguid','rows');


%% load in metadata
load(metaname,'numframes','totalFrameTime','numrepeats','stimONframe','stimOFFframe','laserpower');
numtrials = length(laserpower);
totalFrameTime = totalFrameTime/1000;       %convert into seconds

preStimTime = (stimONframe-1)*totalFrameTime;


[uid,~,stimID_idx] = unique(laserpower);

time = (1:(numframes))*totalFrameTime;
prestimframes = 2:sum(time<preStimTime);
stimframes = stimONframe:stimOFFframe;
iosResponseFrames = round(stimONframe+(.5/totalFrameTime)):round(stimOFFframe+(.5/totalFrameTime));
%% will want to get ONE very good tranformation matrix to map the data from the current session to the master session
% to do this you need to load in some ios data from the current session,
% make an average image; do NOT need to do this if this is the reference
% session:
if ~isReferenceSession
    cd(movdirectory);
    curidx = 1:(max(prestimframes));
    localavios = loadRawStack_v1(rootname,localimagesize,curidx);
    localavios = mean(localavios,3);
    A = localavios;A(~windowmask) = median(A(:));
    B = refim; B(~refim_windowmask) = median(B(:));
    sessionRegister = registerMultipleWidefieldImages_v1_2(A,B);
    localavios_transformed = imwarp(localavios,sessionRegister.Transformation,'OutputView',imref2d([size(refim,1),size(refim,2)]));
    figure; imshowpair(refim,localavios_transformed,"checkerboard");
    clear A B
else
    localavios = avios;
end
%% loop through movies and do the tracking (for .raw formatted image sequences):
cd(movdirectory);
totalframes = numframes*numtrials;
traj = [];
gcampmov_full = [];
parenchymaMask = false(size(bg_ios,1),size(bg_ios,2),numtrials);
roughvesselmask = false(size(bg_ios,1),size(bg_ios,2),numtrials);

avgParenchymaIOS_normalized = zeros(size(parenchymaMask));
for n = 1:numtrials
    curidx = ((n-1)*numframes + 1):(n*numframes);
    iosmov = loadRawStack_v1(rootname,localimagesize,curidx);
    
    iosmov = iosmov-bg_ios;
    if ~isReferenceSession
        [iosmov,curtraj] = trackOneWidefieldMovie_iosonly(iosmov,windowmask,localavios,refim,vesselMasks, sessionRegister,micronsPerPixel);
    else
        [iosmov,curtraj] = trackOneWidefieldMovie_iosonly(iosmov,windowmask,localavios,[],vesselMasks,[],micronsPerPixel);
    end
    [parenchymaMask(:,:,n),roughvesselmask(:,:,n)] = getParenchymaPixels(mean(iosmov,3),windowmask);
    
    for n2 = 1:size(iosmov,3)
    
        temp = iosmov(:,:,n2);
        temp(~parenchymaMask(:,:,n)) = nan;
        temp(~windowmask) = 0;
        iosmov(:,:,n2) = temp;
    end
    temp = mean(iosmov(:,:,prestimframes),3,'omitnan');
    temp = iosmov./repmat(temp,1,1,size(iosmov,3));
    avgParenchymaIOS_normalized(:,:,n) = mean(temp(:,:,iosResponseFrames),3,'omitnan');

    traj = [traj,curtraj];
    disp(['tracked movie # ',num2str(n)])
end
numTrackedSegs = size(traj(1).vesselCent,1);
disp('FINISHED TRACKING DILATIONS')
%% calculate baseline levels for vessel total intensity:
clear iosmov mov curgcampmov curtraj
baselineVesselInten = zeros(numTrackedSegs,1);
for n = 1:length(traj)
    baselineVesselInten = baselineVesselInten+mean(traj(n).vesselTotalInten(:,prestimframes),2,'omitnan');
end
baselineVesselInten = baselineVesselInten/length(traj);

%% calculate  baseline levels for surround total intensity:
baselineSurroundInten = zeros(numTrackedSegs,1);
for n = 1:length(traj)
    baselineSurroundInten = baselineSurroundInten+mean(traj(n).surroundInten(:,prestimframes),2,'omitnan');
end
baselineSurroundInten = baselineSurroundInten/length(traj);

    
%% run processing on gcamp movie AFTER the fact to allow for global baseline correction:
clear iosmov curgcampmov curtraj



%% save data so it's not lost
cd(savedirectory)
if ~isReferenceSession
    save(savename,'baselineVesselInten','traj','parenchymaMask','roughvesselmask','baselineSurroundInten','avgParenchymaIOS_normalized',...
        'laserpower','stimframes','totalFrameTime','sessionRegister');
else
    save(savename,'baselineVesselInten','traj','parenchymaMask','roughvesselmask','baselineSurroundInten','avgParenchymaIOS_normalized',...
        'laserpower','stimframes','totalFrameTime');
end


